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ABSTRACT 

If the favored hierarchical cosmological model is correct, then the Milky Way system should have 
accreted ~ 100 — 200 luminous satellite galaxies in the past ~ 12 Gyr. We model this process using 
a hybrid semi-analytic plus N-body approach which distinguishes explicitly between the evolution of 
light and dark matter in accreted satellites. This distinction is essential to our ability to produce a 
realistic stellar halo, with mass and density profile much like that of our own Galaxy, and a surviving 
satellite population that matches the observed number counts and structural parameter distributions 
of the satellite galaxies of the Milky Way. Our model stellar halos have density profiles which typically 
drop off with radius faster than those of the dark matter. They are assembled from the inside out, with 
the majority of mass (~ 80%) coming from the ~ 15 most massive accretion events. The satellites 
that contribute to the stellar halo have median accretion times of ~ 9 Gyr in the past, while surviving 
satellite systems have median accretion times of ~ 5 Gyr in the past. This implies that stars associated 
with the inner halo should be quite different chemically from stars in surviving satellites and also from 
stars in the outer halo or those liberated in recent disruption events. We briefly discuss the expected 
spatial structure and phase space structure for halos formed in this manner. Searches for this type of 
structure offer a direct test of whether cosmology is indeed hierarchical on small scales. 

Subject headings: Galaxy: evolution — Galaxy: formation — Galaxy:halo — Galaxy: kinematics 
and dynamics — galaxies: dwarf — galaxies: evolution — galaxies: formation 
— galaxies: halos — galaxies: kinematics and dynamics — Local Group — dark 
matter 



1. INTRODUCTION 

There has been a long tradition of searching in the 
stellar halo of our Galaxy for signatures of its forma- 
tion. Stars in the halo provide an important avenue for 
testing theories of galaxy formation because they have 
long orbital time periods, have likely suffered little from 
dissipation effects, and tend to inhabit the outer regions 
of the Galaxy where the potential is relatively smooth 
and slowly evolving. The currently favored Dark Energy 
+ Cold Dark Matter (ACDM) model of structure for- 
mation makes the specific prediction that galaxies like 
the Milky Way form hierarchically, from a series of ac- 
cretion events involving lower-mass systems. This leads 
naturally to the expectation that the stellar halo should 
be formed primarily from disrupted, accreted systems. 
In this work, we develop an explicit, cosmologically- 
motivated model for stellar halo formation using a hy- 
brid N-body plus semi-analytic approach. Set within the 
context of ACDM, we use this model to test the general 
consistency of the hierarchical formation scenario for the 
stellar halo and to provide predictions for upcoming sur- 
veys aimed at probing the accretion history of the Milky 
Way and nearby galaxies^ 

In a classic study, lEggen. Lvnden-Bell. k. Sandagd 
used proper motions and radial velocities of 221 
dwarfs to show that those with lower metallicity (i.e. halo 
stars) tended to move on more highly eccentric orbits. 
They interpreted this trend as a signature of formation 
of the lower metallicity stars during a rapid radial col- 
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lapse. In contrast. iSearle fc ZinnI l)1978j) suggested that 
the wide range of metallicites found in a sample of 19 
globular clusters at a variety of Galactocentric radii in- 
stead indicated that the Galaxy formed from the grad- 
ual agglomeration of many sub-galactic sized pieces. A 
recent analysis of 1203 metal-poor Sola r neighborhood 
stars , selected without kinematic bias (|Chiba fc Beers! 
I200QO, points to the truth being some combination of 
these two pictures: this sample contained a small concen- 
tration of very low metallicity stars on highly eccentric 
orbits (reminiscent of Eggen, Lynden-Bell fc Sandage's, 
1962 work) but otherwise showed no correlation of in- 
creasing orbital eccentricity with decreasing metallicity. 

In the last decade, much more direct evi- 
dence for the lumpy build-up of the Galaxy has 
emerged in the form of clumps of stars in phase- 
space (and, in some cases, metallic ity) bo t h rel- 
atively nearby (Maiewsk i. Munn. fc HawlevI 119961: 
iHelmi. White, de Zeeuw. Z hao 199^ and at much 
larger distances. The most striking example in 
the latter category is the discovery of the Sagit- 
tarius dwarf galaxy flba ta. Gilmore. fc IrwinI Il994t 
Hba ta Gilmore fc^rwi_n_ 19951) — here after Sgr - — and 
its a ssociated trails of debris (see iMaiewski et al. I 
120031 for an overview of the many detections) which 
have now been t r aced entirely around t he Galaxy 
pbata et all 120011: IMaiewski et al. I 120031) . Large 
scale surve ys of the s tellar halo a re now under- 
way (Maiewski. Osthcimcr. Kunke^^^Pattgggn| 
20iB iMorrison et al. .2000: Yannv et alJ I200nt 
Ivezic et all l2000t iNewberg et all 120021) . and have 
uncovered additional structures, not associated 
with Sgr ( Ncwberg et al. 2003; Martin et al. 200l 
iRocha-Pinto et al.l l2004|) . Moreover, recent advances 
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in instrumentation are now permitting searches for 
and discoveries of analogous structures around other 
galaxies in th e form of overdensities in integrated light 
fehang et al.llT9 98: Zheng ct al. 1999; Forbes ct al. 2003; 
iPohlen et al II2003I) or, in the case of M31, star counts 
ilbata et all 120011: lE^guson et alJ 120021: iZucker et all 
|2004() . Given this plethora of discoveries, there can 
be little doubt that the accretion of satellites has 
been an important contributor to the formation of our 
and other stellar halos . In add ition, both theo retical 
iAbadi et a l. 200 3 : Bro ok et a £ 2GG5a; Rober tson et al. 
222^ and observational (Gilmore, Wvs^^^^Ngmj 
2f)MlYa, nnv et a,l. H2003 : Tbata, et a,l.ll2003t iGrane et ab 
2003; Ro cha-Pinto. Maie wski. Skrutskie. fc Cranell200a 
Frinchabov et all l2004l:THelm'i et all l200g ) work is be- 
ginning to suggest that some significant fraction of the 
Galactic disk could also have been formed this way. 

All of the above discoveries are in qualitative agree- 
ment with the expectation s of hierarchical struc- 
ture formation fPeeblesI 119651 IPress fc SchechteJ [1971 
[Blumenthal et al, 1984). As the prevailing variant of this 
picture, ACDM is remarkably successful at reproducing a 
wide ran ge of observ ations, especiall y on large scal es (e.g. 
pisenstei n et al.ll200 5: Mail er et alJl20 05: Tcgmar k et al.l 
l2004HSpereel et alJi2003l: iPercival et al.ii2002^) . On sub- 
galactic scales, however, the agreement b etween theory 
and observation is not as obvious (e.g. Simon et al. 200"! 
iKazantzidis et alJl2004 iD'Onghia fc Burkert..20041 . In- 
deed, the problems explaining galaxy rotation curve 
data, dwarf galaxy counts, and galaxy disk sizes have 
lead some to suggest modifications to the standard 
paradig m, including an allowance for warm dark mat- 
ter (e.g. lSommer -Larsen et al.l2004D . early-decaying dark 
matter lIKa plinghat 2005 ), or non-standard inflation 
ijZentnerfc BuUockl 120021 12003(1 . These modifications 
generally suppress fluctuation amplitudes on small scales, 
driving sub-galactic structure formation towards a more 
monolithic, non-hierarchical collapse. These issues bring 
to sharper focus a fundamental question in cosmology 
today: is structure formation truly hierarchical on small 
scales? Stellar halo surveys offer powerful data sets for 
directly answering this question. 

Numerical simulations of individual satellites dis- 
rupting about parent galaxies can in many cases 
provide convincing similarities to the observed phase- 
space lumps. These models allow the observations 
to be interpreted in terms of t he mass and orbit 
of the pro genitor satellite (e.g.. iVelazaiiez fc Whitel 
19951: I Johnston. Soergel. fc HernauistI 119951 

^ghngim^SiEurdssom'&Hem^^ et al.' 

1999t iHel mi fc White 2001; Helmi et al. 2003; 
Law. Johnston fc Maiewski! |2004|) . and even the 
potenti al of the galaxy in which it is orbit- 

^ HJohnston. Zhao. Snergel. fc HernnuisH [l^ 

Miirali fc D ubinsk il [l99a llbata et alJ I2001I l2004t 
Johnston. L aw fc Maiewskill2005|) . Nevertheless, a true 
test of hierarchical galaxy formation will require robust 
predictions for the frequency and character of the 
expected phase space structure of the halo. 

Going beyond qualitative statements to model the 
full stellar halo (including substructure) within a cos- 
mological context is non-trivial. The largest contribu- 
tor of substructure to our own halo is Sgr, estimated 
to have a currently-bound mass of order 3 x IO^Mq 



l|Law. Johnston fc Maiewskil |2004|) . Even the highest 
resolution cosmological N-body simulations would not 
resolve such an object with more than a few hundred 
particles, which would permit only a poor represen- 
tation of the phase-spa ce structure of its debris (see 
IHelmi. White, fc Snringel 2003. for an example of what 
can currently be done in this field). Such simulations 
are computationally intensive, so the cost of examining 
more than a handful of halos is prohibitive and it is diffi- 
cult to make statements about the variance of properties 
of halos that might be seen in a large sample of galax- 
ies. Moreover, such simulations in general only follow 
the dark matter component of each galaxy not the stellar 
component. In their studies of thick disk and inner halo 
formation. Brook and collaborators ifB rook et al.ll200'^ 
l2004a b. 2005a b) have modeled the stellar components 
directly by simulating the evolution of individual galaxies 
as isolated spheres of dark matter and gas with small- 
scale density fluctuations superimposed to account for 
the large-scale cosmology. However, their sample size 
remains small and, though they are able to make gen- 
eral statements about the properties of their stellar ha- 
los, their resolution would prohibit a detailed phase-space 
analysis. 

An alternative is to take an analytic 
or semi-analytic appr oach to halo build- 
ing (e.g. [B ullock. Kra vtsoy. fc W einberg 200^ 
I Johnston. Sackett. fc BuUockl l200l[ .Tavlor. .200^ . 
This allows the production of many halos, and the 
potential of including prescriptions to follow the stars 
separately from the dark matter. However, such 
techniques use only approximate descriptions of the 
dynamics and are unable to follow the fine details of the 
phase-space structure accurately. 

In this study we develop a hybrid scheme, which draws 
on the strengths of each of the former techniques to 
build high resolution, full phase-space models of a sta- 
tistical sample of stellar halos. Our approach is to vastly 
decrease the computational cost of a full cosmological 
simulation by modeling only those accretion events that 
contribute directly to the stellar halo in detail with N- 
body simulations, and to represent the rest of the galaxy 
with smoothly-evolving analytic functions. The baryonic 
component of each contributing event is followed using 
semi-analytic prescriptions. 

The purpose of this paper is to describe our method, its 
strengths and limitations (§2), present the results of tests 
of the consistency of our models with general properties 
of galaxies and their satellite systems (§3) and outline 
some implications (§4). We summarize the conclusions 
in §5. In further work we will go on to compare the 
full phase-space structure of our halos in detail to obser- 
vations and to examine the evolution of dark and light 
matter in satellite galaxies after their accretion. 

2. METHODS 

Our methods can be broadly separated into: (I) a sim- 
ulation phase, which follows the phase-space evolution 
of the dark matter; and (II) a prescription phase, which 
embeds a stellar mass with each dark matter particle. 
Specifically: 

Phase I: Simulations 
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A We generate merger trees for our parent galaxies 
using the method outlined in iSomerville fc Kolatd 
(|1999l) based o n the Ex tended-Press-Schechter 
(EPS) formalism l|Lacev fc Cole 1991 — see §2.1)). 

B For each event in step I A, we run a high-resolution 
N-body simulation that tracks the evolution of the 
dark matter component of a satellite disrupting 
within an analytic, time-dependent, parent galaxy 
-I- host halo potential (see §2.2). 

Phase II: Prescriptions 

A We follow the gas accretion history of each 
satellite prior to falling into the parent and 
track its star-formation rate using cosmologically- 
motivated, semi-analytic prescriptions (see §2.3). 

B We embed the stellar components generated in step 
IIA within each dark matter satellite by assigning 
a variable mass-to-light ratio to every particle that 
is tracked in the (Phase I) N-body simulations (see 
§2-4). 

We consider the two-phase approach a necessary and 
acceptable simplification since it allows us to separate 
well-understood and justified approximations in Phase 
I from prescriptions that can be adjusted and refined 
during Phase II. In addition, this separation allows us 
to save computational time and use just one set of dark 
matter simulations to explore the effect of varying the 
details of how baryons are assigned to each satellite. A 
more complete discussion of the strengths and limitations 
of our scheme is given in §2.5. 

2.1. Cosmological Framework 

Throughout this work we assume a ACDM cosmology 
with = 0.3, f^A = 0.7, f^b/i^ = 0.024, h = 0.7, and 
(Tg = 0.9. The implied baryon fraction is rif,/r2ni = 0.16. 

We focus on the formation of stellar halos for "Milky- 
Way" type galaxies. In all cases our z ~ host dark 
matter halos have virial masses Mvi,. = 1.4 X lO^^^fo, 
corresponding virial radii i?vir,o = 282 kpc, and virial 
velocities V^ir = 144km s~^. The quantities M^-„ and 
i?vir are related by 

47r 

Mvir = —pM[z)Kh{z)Rl;r, (1) 

where pm is the average matter density of the Universe 
and Avir is the "virial overdensity" . In the cosmology 
considered here, Avir(-z = 0) ~ 337, and Avir 178 
at 1 (jBrvan fc Norma'rilll998() . The virial velocity is 

defined as Kir = \/GMvir/i?vir- 

We generate a total of eleven random realizations of 
stellar halos. General properties of all eleven are sum- 
marized in Table Q Any variations in our results for 
stellar halos among these are determined by differences 
in their accretion histories. In all subsequent figures we 
present results for four stellar halos (1,2,6, and 9) chosen 
to span the range of properties seen in our full sample. 

2.1.1. Semi- analytic accretion histories 

We track the mass accretion and satellite acqui- 
sition of each parent galaxy by constructing merger 



trees using the statistica l Monte Carlo method of 
ISomerville fc KolatU il999D based on the EPS formal- 
ism (ILacev fc Colelll993D . This method gives us a record 
of the the masses and accretion times of all satellite 
halos and hence allows us to follow the mass accre- 
tion history of each parent as a function of lookback 
time. We explicitly note all satellites more massive 
than Mmin = 5 x IQ^Mq and treat all smaller accre- 
tion events as diffuse mass accretion. Column 2 of Table 
^ lists the total number of such events for e ach simu- 
lated halo. For furth er deta i ls see | Lacev fc C ole (1993); 
ISomerville fc KolatJ l(T999|) : IZentner fc Mlock (2003).' 
Four examples of the cumulative mass accretion histories 
of parent galaxies generated in this manner are shown by 
the (jagged) solid lines in Figure ^ 

2.1.2. Satellite orbits 

Upon accretion onto the host, each satellite is assigned 
an initial orbital energy based on the range of b inding en- 
ergies observed in cosmological simulations (|Klvpin et alJ 
I1999|) . This is done by placing each satellite on an ini- 
tial orbit of energy equal to the energy of a circular orbit 
of radius i?circ = V^viv, with 77 drawn randomly from a 
uniform distribution on the interval [0.4,0.8]. Here i?vir 
is the virial radius of the host halo at the time of accre- 
tion. We assign each subhalo an initial specific angular 
momentum J = eJcirc, where Jcirc is the specific angu- 
lar momentum of the aforementioned circular orbit and 
e is the orbital circularity, which takes a value between 
and 1. We ch oose e from the binned dis tribution shown 
in Figure 2 of iZentner fc Bullockl l(200l . which was de- 
signed to match the c osmol ogical N-body simulation re- 
sults of 'Ghigna et aT (1998*), and is similar to the circu- 
larity di stributions found in more recent N-body analyses 
l|Zentner et al.1120041: IBenson 2005). Finally, the plane of 
the orbit is drawn from a uniform distribution covering 
the halo sphere. 

2.1.3. Dark matter density distributions 

We model all satellite and pa rent halos with the spheri- 
cally averaged density profile of lNavarro. Frenk fc Whit^ 
(,1996) (NEW): 

where rhaio (= in NEW) is the characteristic inner 
scale radius of the halo. The normalization, ps, is set 
by the requirement that the mass interior to i?vir be 
equal to Mvir. The value of rhaio is usually character- 
ized in terms of of the halo "concentration" parameter: 
c = Rvir/ rhaio- The implied maximum circular veloc- 
ity for this profile occurs at a radius rmax — 2.15rhaio 
and takes the value Knax — 0.466Vvir-F(c), where F{c) — 
v/c/[ln(l + c)-c/(l + c)]. 

For satell ites, we set the valu e of c using the simulation 
results of lBuUock et al.l l|200lD and the corresponding re- 
lationship between halo mass, redshift, and concentra- 
tion summarized by their analytic model. The median 
c relation for halos of mass Mvir at redshift z is given 
approximately by 
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although in practice we use t he full analytic model dis- 
cussed in lBullock eFall l|2n0lD . 

For parent halos, we allow their concentrations 
to evolve self-consistently as their virial masses in- 
crease, as has been s een in the N-body simulations of 
iWechsler et alJ l)2fl02D . Rather than represent the halo 
growth as a series of discrete accretion events, we smooth 
over the Monte Carlo EPS merger tree by fitting the fol- 
lowing functional form to the Monte Carlo mass accretion 
history for each halo: 



Mvir(a) = Afvir(ao) cxp 



-2a. - 1 
V a 



(4) 



Here a = (1 + z)^^ is the expansion factor, and a. is 
the fitting parameter, corresponding to the value of the 
expansion factor at a characteristic "epoch of collapse" . 
IWechsler et alJ 1)20021) demonstrated that the value of Oc 
connects in a one-to-one fa shion with t he halo concen- 
tration parameter t Wcchsle r et al.ll2002|) : 



c(a) = 5.1 — . 



(5) 



Halos that form earlier (smaller a.'s) are more concen- 
trated. 

Example fits to four of our halo mass accretion histories 
are shown by the smooth solid lines in Figure 1. The a. 
values for each of the halos in this analysis are listed in 
the third column of Table 1. Typical host halos in our 
sample have c ~ 14 at z = 0, scale radii rhaio — 20kpc, 
and maximum circular velocities l^nax — 190km s~^. 

2.2. N-body simulations of dark matter evolution 

Having determined the mass, accretion time and orbit 
of each satellite (§2.1.1 and §2.1.2), and the evolution 
the potential into which it is falling (§2.1.3), we next 
run individual N-body simulations to track the dynami- 
cal evolution of each satellite halo separately. We follow 
only those that contain a significant stellar component 
(see §2.3 below). In practice, this restricts our analysis 
to satellite halos more massive than Afvir^ IO^Mq — 
the number of such satellites infalling into each parent 
is listed in column 5 of Table Based on our star- 
formation prescription discussed in §2.3, systems smaller 
than this never contain an appreciable number of stars 
and thus don't contribute significantly to the stellar halo. 

2.2.1. The parent galaxy potential 

The parent galaxy is represented by a three-component 
bulge/disk/dark halo potential which we allow to evolve 
with time as the halo accretes mass. The (spherically- 
symmetric) dark halo potential at each epoch a is given 
by the NFW potential generated by the dark matter dis- 
tribution in equation 



$halo(T') 



1 
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T-halo {r/rha.lo) \ ''halo 



1 



(6) 



equation where Mhaio = A^haio(a) and rhaio = ?'haio(a) 
are the instantaneous mass and length scales of the halo 
respectively. The halo mass scale is related to the virial 
mass via 



ln(c-|-l)-c/(c+l)' 



(7) 



The disk and bulge are assumed to grow in mass and 
scale with the halo virial mass and radius: 



disk 



(i?, z) = - 



^disk 



sphere 



(r) 
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(9) 



where Afdisk(a) = 1-0 x lO"(Afvir/A/vir,o)M0, 

Arsphorc(a) = 3.4 X lOlO(A/vir/Afvir,o)M0, iJdisk = 

6.5(rvir/rvir,o) kpc, Zdisk = 0.26(rvir/rvir,o) kpc and 

^sphere = 0.7(rvir/rvir,o) kpC. 

2.2.2. Satellite initial conditions 

We use 10^ particles to represent the dark matter in 
each accreted satellite. Particles are initially distributed 
as an isotropic NFW model, with mass and scale cho- 
sen as described in §2.1.2. The phase-space distribution 
function is derived by integrating over the density and 
potential distributions 



with p 



d'^p d^ 



1 



d^^ Ve^* 



dp 



(10) 



and where \E' = — <&nfw + '^'o is the rel- 
ative potential (such that * ^ as r ^ oo) and e = 
— /2 is the relative energy fsee lBinnev Sz Tremaind 
IT987i for discussion). This distribution function is 
used (in tabulated form) to generate a random realiza- 
tion. This ensures a stable satellite configuration — 
initial conditions generated by instead assuming a lo- 
cal Maxwellian velocity distribution have been shown to 
evolve (Kazantzidis. Magorrian. & Moore 2004). Given 
/(e), the differential energy distribution follows in a 
straightforward manner from the density of states, g{e), 



/(e)5(e), 5(e) = 1671^ / y^2{¥^ydr, 



dM 

~dr 

(11) 

where is the largest energy that can be reached by a 
star of relative energy e. The differential energy distribu- 
tion for our initial halo is shown by the solid histogram in 
Figure 13 We see that the majority of the (dark matter) 
material in an infalling satellite is quite loosely bound. 

Rather than generating a unique /(e) and particle dis- 
tribution for each satellite in each accretion history, a 
single initial conditions file with unit mass and scale, and 
outer radius i?out = 35rhaio {— 35 in our units) is used 
for all simulations with masses and scales appropriately 
rescaled for each run. Since all of our accreted satellites 
have concentrations c < 35, our set up effectively allows 
each accreted satellite's mass profile to extend beyond 
its virial radius for several scale lengths. We do not ex- 
pect this simplification to significantly affect our results 
because the the light matter is always embedded at the 
very central regions of the halo (r*^ rhaio) and the outer 
material is always quickly stripped away from the outer 
parts of the halos upon accretion. 

In §2.4 we discuss our method of "embedding" star 
particles within the cores the accreted satellite dark ha- 
los. 
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2.2.3. Satellite evolution 

The mutual interactions of the satelhte particles 
are calculated using a bas is function expansion code 
ijHernauist fc OstrikedlTOQ^ . The initial conditions file 
for the satellite is allowed to relax in isolation for ten 
dynamical times using this code to confirm stability. For 
each accretion event a single simulation is run, following 
the evolution of the relaxed satellite under the influence 
of its own and the parent galaxy's potential, for the time 
since it was accreted (as generated by methods in §2.1.1) 
along the orbit chosen at random from the distribution 
discussed in §2.1.2. (Note that simulations of satellite 
accretions in static NFW potentials using thi s code pro- 
duced results identical to those reported in .Havashi et al.l 

l2nnl 

Using this approach, the satellites are not influenced 
by each other, other than through the smooth growth 
of the parent galaxy potential. Nor does the par- 
ent galaxy react to the satellite directly. In order to 
mimic the expected decay of the satellite orbits due to 
dynamical friction (i.e. the interaction with the par- 
ent), we include a drag term on all particles within 
two tidal ra dii rt.idp of the satellite's center, of the form 
proposed by iHashimoto. Funa to. fc Makind ||2003D an d 
modified for NFW hosts bv Zcntncr fc BullocU l)2003j) . 
This approach includes a slight modification to the stan- 
dard Chandrasekhar d ynamical friction formula (e.g. 
iBinnev fc Trem aine"1987) . The tidal radius rtido is calcu- 
lated from the instantaneous bound mass of the satellite 
TTisat, the distance r of the satellite to the center of the 
parent galaxy and the mass of the parent galaxy within 
that radius as rtidc — ?'(™sat/Afr)^^'^- 

2.2.4. Increasing phase-space resolution with test 
particles 

In this study, we are most interested in following the 
phase-space evolution of the stellar material associated 
with each satellite. This is assumed to be embedded 
deep within each dark matter halo (see §2.4) — typically 
only of order 10'* of the N-body particles in each satellite 
have any light associated with them at all. In order to 
increase the statistical accuracy our analysis we sample 
the inner 12% of the energy distribution with an addi- 
tional 1.2 X 10^ test particles. This does not increase the 
dynamic range our simulation, but does allow us to more 
finely resolve the low surface brightness features we are 
interested in with only a modest increase in computa- 
tional cost: we gain a factor of 10 in particle resolution 
with an increase of ~25% in computing time. In this 
paper, we have used test particles only in generating the 
images shown in Figures E| - 1161 

2.3. Following the satellites ' baryonic component 

We follow each satellite's baryonic component using 
the expected mass accretion history of each satellite halo 
(prior to falling into the parent galaxy) in order to track 
the inflow of gas. The gas mass is then used to de- 
termine the instantaneous star formation rate and to 
track the buildup of stars within each halo. The physics 
of galaxy formation is poorly understood, and any at- 
tempt to model star formation and gas inflow into galax- 
ies (whether semi-analytic or hydrodynamic) necessarily 
require free parameters. Our own prescription requires 



three "free" parameters: Zro, the redshift of reionization 
(see H2.3.1|l : /gas, the fraction of baryonic material in the 
form of cold gas (i.e. capable of forming stars) that re- 
mains bound to each satellite at accretion (see H2.3.2|l : 
and ti,, the globally-averaged star formation timescale 
(see ^TTT^ . 

In the following subsections we describe how these pa- 
rameters enter into our prescriptions, and choose a value 
of /gas consistent with observations. In §3 we go on to 
demonstrate that the observed characteristics of the stel- 
lar halo (e.g. its mass, and radial profile) and the Milky 
Way's satellite system (e.g. their number and distribu- 
tion in structural parameters) provide strong constraints 
on the remaining free parameters and hence the efficiency 
of star formation in low-mass dark matter halos in gen- 
eral. 

2.3.1. Reionization 

Any attempt to model stellar halo buildup within the 
context of ACDM must first confront the so-called "miss- 
ing satellite problem" — the apparent over-prediction of 
low-mass halos compared to the abundance of satellite 
galaxies around the Milky Way and M31. For example, 
there are eleven know satellites of the Milky Way — nine 
classified as dwarf spheroidal and two as dwarf irregulars 
— yet numerical work predicts several hundred dark mat- 
ter satellite halos in a s imilar mass range ((Klvpin et alJ 
ll999HMoore et alJl999D . It is quite likely that our inven- 
tory of stellar satellites is not complete given the lumi- 
nosity and surface brightness limits of prior searches (as 
the recent discove ry of the Ursa M inor dwarf spheroidal 
demonstrates, see lWillman et aD[2005l) . but incomplete- 
ness is no t seen as a viable solu tion for a problem of this 
scale Csee lWillman et al]l2004i for a discussion). 

The simplest solution to this problem is to postulate 
that only a small fraction of the satellite halos orbiting 
the Milky Way host an observable galaxy. In this work, 
we solve the missing satellite problem using t he sug- 
gestion of lBullock. Kravtsov. fc Weinberg 1 1)2000() . which 
maintains that only the ~ 10% of low-mass galaxies 
(Knax < 30km s^*) that had accreted a substantial frac- 
tion of their gas before th e epoch o f reionization host 
obser vable galaxies (see alsolChiu et ah 2001: Somcrvill^ 
l200llBenson et al.ll200l iKravtsov et alJl2004|) . The key 
assumption is that after the redshift of hydrogen reion- 
ization, Zie, gas accretion is suppressed in halos with 
Vmax < 50km s~*, and completely stopped in halos with 
Knax < 30km s^*. These th r eshold s foll ow from t he re- 
sults of iThon] fc Weinberg! llT993) and ICnedinI Eml 
who used hydrodynamic simulations to show that gas 
accretion in low-mass halos is indeed suppressed in the 
presence of an ionizing background. 

We also impose a low-mass cutoff for tracking galaxy 
formation in satellite halos with VJnax < 15km s~*. Two 
processes and one practical consideration motivate us 
to ignore galaxy formation in these tiny halos: first, 
photo-evaporation acts to eliminate any gas that was ac- 
creted before reion ization in halos with Vpiax^ 15km s~* 
l|Barkana fc Loebl ll999: Sh aviv fc Dekell2n03D : second, 
the cooling barrier below virial temperatures of ^ lO^K 
(corresponding to Knax ^ 16kms~^) prevents any gas 
that could rem ain bound to t hese halos from cooling and 
forming stars ijKepner et al.i,1997; ..Dekel fc Woo 2003) : 
finally, even if we were to allow star formation in these 
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systems, their contribution to the stellar halo mass would 
be negligible. Once we are more confident of our inven- 
tory of the lowest luminosity and lowest surface bright- 
ness satellites ( Willman et al. 2003)of the Milky Way we 
should be able to confirm these physical arguments with 
observational constraints. 

The epoch of reionization Zre determines the numbers 
of galaxies that have collapsed in each of the above Vmax 
limits, and hence the number of luminous satellites that 
will be accreted, whether they disrupt to form the stellar 
halo or survive to form the Galaxy's satellite system. We 
discuss limits on this parameter in §3.1.1. 

2.3.2. Gas accretion following reionization 

The virial mass of each satellite, M^f^, at the time of 
its accretion, Cac, is set by our merger tree initial con- 
ditions (§2.1.1). We assume that each satellite halo has 
had a mass accumulation history set by Equation0]up to 
the time of its merger into the "Milky Way" host, with 
o-a = flac- After accretion, all mass accumulation onto the 
satellite is truncated (see §2.3.3). For massive satellites, 
Knax > 50km s~^, we set Oc in Equation^using the satel- 
lite's mass-defined concentration parameter via E quation 
n fsee IRu]1ock7nni2nnTl IWechsler et a,1 .11200^. This 
provides a "typical" formation history for each satellite. 
For low-mass satellites, we are necessarily interested in 
where Oc falls in the distribution of halo formation epochs 
because this determines the fraction of mass in place at 
reionization. Therefore, if Vmn^ < 50km s~^, we use the 



methods of iLacev fc Cold l[l993|) in order to derive the 
fraction of the satellite's mass that was in place at the 
epoch of reionization, z^e, and use this to set the value 
of ttc- Given Qc for each satellite, we determine the in- 
stantaneous accretion rate of dark matter h{t) in to this 
system as a function of cosmic time via 



d M^?* 

h(t) = 
^ ' dt 



(12) 



In the absence of radiative feedback effects, cooling is 
extremely efficien t in pre-merged satellite s of the size we 
consider (see, e.g. lMaller fc Bullockll2004|l . Therefore we 
expect the cold gas inflow rate to track the dark matter 
accretion rate, h{t) — at least in the absence of the effects 
of reionization — and take it to be 

C/gasMi-^in). (13) 

The time lag within h{t) accounts for the finite time it 
takes for gas to settle into the center of the satellite af- 
ter being accreted. We assume this occurs in roughly a 
halo orbital time at the virial radius: tin = T^Rh/Vvii — 
6 Gyr (1 + z)~^/^. We have introduced the constant C 
in order to account for the suppression of gas accretion 
in low- mass halos (as alluded to in tl2.3.1|l . Before the 
epoch of reionization, we set C 1 for all galaxies. For 
systems with Vniax > 50km s^^, C = 1 at all times. After 
reionization, C = in systems with Vmax < 30km s~^, 
and C varies linearly in Vmax between and 1 if Vm 



falls between 30 and 50km s ^ (see iThoul fc Weinberd 

The fraction of mass in each satellite in the form of 
cold, accreting baryons, /gas, determines the total stellar 
mass plus cold gas mass associated with each dark matter 
halo. In what follows, we adopt /gas ~ 0.02, which is an 
upper limit on the range of cold baryonic mass fraction 
in observed galaxies (Bell et al. 2003). 



2.3.3. Star Formation 

If we assume that cold gas forms stars over a timescale 
t^,, then the evolution of stellar mass M^, and cold gas 
mass Mgas follows a simple set of equations: 



dA4 M, 
"dT 

dA/gas 



gas 



dAf. 



+ C/gas/l(t-<in). 



(14) 

dt dt ' (15) 

For simplicity, the star formation is truncated soon 
after each satellite halo is accreted onto the Milky 
Way host. Physically, this could result from 
gas loss via ram-p ressure stripping fr om the back- 
groun d hot gas halo iLi n fc F ^ bcr 1983; Moore & Davr 
isha^ l2r 



19941 iB litz & Elobisha^ |22 12 jailer & Bullock 200 
Maver et al.«2005^1 . This model is broadly consistent with 
observations that demonstrate that the gas fraction in 
satellites of the Milky Way and Andromeda is typically 
far less than that in field dwarfs in the Local Group, 
as illustrated by the separation of the open (satellites) 
and filled (fiel d dwarfs) symbols in Figure El (plotting 
data taken from lGrebel. Gallagher, fc* Harbcck 2003). Of 
course, this assumption is over-simplified, but it allows 
us to capture in general both the expectations of the hi- 
erarchical picture and the observational constraints. We 
note that this is likely a bad approximation for massive 
satellites, whose deep potential wells will tend to resist 
the effects of ram-pressure stripping. However, we ex- 
pect that this will have little impact on our stellar halo 
predictions, since most of the stellar halo is formed from 
satellites that are accreted early and destroyed soon af- 
ter. 

The star formation timescale, determines the star to 
cold gas fraction in each satellite upon accretion and, for 
a given value of /gas, total stellar luminosity associated 
with each surviving satellite and the stellar halo. We 
discuss limits on this parameter in §3.1.2. 

2.4. Embedding baryons within the dark matter satellites 

We model the evolution of a two-component popula- 
tion of stellar matter and dark matter in each satellite 
by associating stellar matter with the more tightly bound 
material in the halo. As discussed in i l2.1.3l the mass 
profile of the satellite is assumed to take the NFW form. 
Mass-to-light ratios for each particle are picked based on 
the particle energy in order to produce a realistic stellar 
profile for a dwarf galaxy. 

A phenomenologically-motivated approximation for 
the stellar distribution in dwarf galaxies is the spheri- 
cally symmetric King profile (King 1962): 



K cos-Hx) r -\ _ 1 + (r/r,)^ 

J' "-n-(nAe)2- 

(16) 

The core radius is Tc and rt is the tidal radius, where 
p^(r > rt) — 0. The normalization, K, is set by the 
average density of the satellite, determined by its mass 
H2.3.3(l and size scales (discussed below). 

For each satellite, we assume a stellar mass to light 
ratio of Af*/_Lv = 2, and use the stellar mass calculated 
m in order to assign a median King core radius 

/ L \°-^^ 
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where throughout L^, is assumed to be the V-band stehar 
luminosity. We allow scatter about the relation using a 
uniform logarithmic deviate between —0.3 < A log^Q L < 
0.3. This slope and normalization was determined by 
least-square fit to the luminosity and core size corre- 
lation for the dwarf satellite data presented in Mateo 
(1998), and the scatter was determined by a "by-eye" 
comparison to the scatter in the data about the relation. 
Our adopted relation between Tj, and is also consistent 
with the relevant projection of the fundamental plane for 
dwarf galaxies (e.g. Kormendy 1985). For all satellites 
we adopt n/rc = 10. 

Assuming isotropic orbits for the stars and that the 
gravitational potential is completely dominated by the 
dark matter, the stellar energy distribution function cor- 
responding to the King profile /*(e) is determined by 
setting p — Pi, and 'J = — $nfw + in equation ((11])) . 
The mass-to-light ratio of a particle of energy e is then 
simply /^(e)//(e) = (dM/de)^/(dM/de). Three exam- 
ples are given Figure El 

2.5. Limitations of our method 

The main limitation of our method is that it only 
follows the smooth growth of the parent potential an- 
alytically — the satellite/satellite interactions and reac- 
tion of the parent to the satellite are not modeled self- 
consistently. Hence we do not anticipate following the 
evolution of the field or satellite particles during a major 
or even minor merger event with great accuracy. Given 
this limitation, we only simulate the accretion histories of 
halos generated from the Monte Carlo merger tree code 
that have not suffered a significant merger (>10 % of 
the parent halo mass) in the recent past (<7 Gyr) — 
11 of the 20 accretion histories generated met this crite- 
rion. In addition, we consider results from simulations 
of accretion events that have occurred prior to the last 
significant merger to be less reliable. We label the halos 
used in this work 1-11. The five left-hand columns of 
Table n summarize the properties of the simulations run 
for each halo. 

Even with these restrictions, we consider our study to 
be a useful approach for exploring substructure in galaxy 
halos because: (i) the highest surface brightness features 
in halos are likely to have come from recent events, whose 
debris has had a shorter time to phase-mix and/or be 
dispersed by oscillations in the parent galaxy potential; 
and (ii) substructure should be more readily detectable 
around spiral (rather than elliptical) galaxies because 
their stellar distributions are less extended — the exis- 
tence of disks in spirals suggests that these are the ones 
the more quiescent accretion histories. 

3. RESULTS I: TESTS OF THE MODEL 

As outlined in §2.5, our method most accurately fol- 
lows the phase-space evolution of debris from accretion 
events that occur during relatively quiescent times in a 
galaxy's history (which we define as being after the last 
>10% merger event). In future work we concentrate on 
those events. In this paper, we analyze the results from 
simulations of the full accretion histories of our halos. 
While not accurate in following the phase-space proper- 
ties of debris material from events occurring before the 
epoch of major merging, the fact that these systems are 
disrupted is predicted robustly, and we are able to record 



the time of disruption and the cumulative mass in those 
disrupted events as well. 

In what follows, we first constrain the remaining free 
parameters z^o (§3.1.1) and ti, (§3.1.2) (with /gas = 0.02), 
by requiring that the general properties of our surviv- 
ing satellite populations are consistent with those of the 
Milky Way's own satellites. We then go on to demon- 
strate that these parameter choices naturally produce the 
observed distributions in and correlations of the struc- 
tural parameters of surviving satellites (see §3.2.1 and 
3.2.2), as well as stellar halos with total luminosity and 
radial profiles consistent with the Milky Way (§3.2.3). 

3.L Primary constraints on parameters 
3.1.1. Satellite number counts 

As described in §2.3.1, we have chosen to solve the 
missing satellite problem by suppressing gas accretion 
in small halos after the epoch of reionization, Zio, and 
suppressing gas accumulation all together in satellites 
smaller than 15km s~^. The number of satellites that 
host stars is then set by choosing z^e. In the work pre- 
sented in this paper, we assume that reionization oc- 
curred at a redshift Zj-c — 10 or at a lookback time of 
13Gyr. The fifth column of Table Ogives the number of 
luminous satellites accreted over the lifetime of each halo 
and the sixth column gives the number of luminous satel- 
lites that survive disruption in each. (The numbers in 
brackets are for those events since the last >10% merger.) 
We see that our reionization prescription leads to agree- 
ment within a factor of ~ 2 with the number of satellites 
observed orbiting the Milky Way. Our results are roughly 
insensitive to this choice as long as 8^ Zic^ 15. 

3.1.2. Infalling satellite gas content 

When reviewing the properties of Local Group 
dwarf galaxies it is striking that — with the no- 
table exceptions of the Large and Small Magellanic 
Clouds (hereafter LMC and SMC) — satellites of 
the Milky Way and Andromeda galaxies are exceed- 
ingly gas-poor compared to their field c ounterparts 
(■Mateo 1998: " Grebel. Gallagher. & Harbeckll2003D . Fig- 
ure 121 emphasizes this point by plotting the V-band 
luminosity vs gas fraction f r om t he compilation by 
iGrebel. Gallagher, fc Harbeckl ()2003D for satellites (open 
squares) and field dwarfs (filled squares). We see that 
field dwarfs tend to have Mhi/Lv — 0.3 — 3, whereas 
satellite dwarfs have gas fractions ~ 0.001 — 0.1. 

While our star formation model assumes that most of 
the gas in accreted dwarfs is lost shortly after a dwarfs be- 
comes a satellite galaxy, consistency with the field dwarf 
population requires that the most recent events in our 
simulations have gas-to-star ratios of order unity imme- 
diately prior to their accretion. This requirement forces 
us to choose a long star formation timescale, = 15 Gyr, 
comparable to the Hubble time. Figure^shows the ratio 
A/gas /^v each satellite at the time it was accreted for our 
four example halos. The clear trend with accretion time 
follows because early accreted systems have not had time 
to turn their gas into stars. Solid points indicate satel- 
lites that survive until the present day. We see that the 
most recently accreted systems (taccr ^1 — 2 Gyr, those 
that should correspond most closely with true "field" 
dwarfs today) have Mgas/i* ^1 — 2, which is in reason- 
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able agreement with the gas content of field dwarfs. The 
points along the lower edge of the trend have lower gas 
fractions at a fixed accretion time because they stopped 
accreting gas at reionization (see §2.3.1). 

Our choice of = 15Gyr is much longer than is typ- 
ical for semi-analytic prescriptio ns of galaxy formation 
set w ithin the CDM context fe.g.. lSomerville fc PrimacS 
Il999|) . but these usually focus on much larger galaxies 
than the dwarfs we focus on here, where star formation 
is likely to have proceeded more efficiently. Observations 
suggest that the dwarf spheroidal satellites of the Milky 
Way have rather bursty, sporadic star formatio n histo- 
ries, with rec ent star fo rmation in some cases 
pObO; Smcckcr -Hane fc Mc William 1999; Gallar t et all 
^99). This effectively demands that the star formation 
timescales must be long in these systems: our model can 
be viewed as smoothing over these histories with an av- 
erage low-level of star formation. 

Note that we do not explicitly include supernova feed- 
back in our star formation histories, but it is implicitly 
included by requiring a very low level of efficiency in our 
model (i.e. a large value of t^^). In two companion papers 
we do include the effects of feedback (accounting for both 
gas gained due to mass loss from stars during normal 
stellar evolutionary phases and gas lost via winds driven 
by supernovae) in order to accurately model chemical 
enrichment in our ac creted satellites (jRobertson et al.l 
l2005HFont et alJ2005|) . With feedback included, a choice 
of ti, = 6.75Gyr provides nearly identical distributions of 
gas and stellar mass in satellite galaxies as does our non- 
feedback choice of ~ 15Gyr. 

3.2. Verification of Model's Validity 

3.2.1. Distributions in satellite structural parameters 

Figure [S] shows histograms of the fractional number of 
satellites as a function of central surface brightness /xg, 
total luminosity and central line-of-sight velocity dis- 
persion (Tj, for the Milky Way dwarf spheroidal satellites 
in solid lines. The dashed lines represent our simulated 
distribution of surviving satellite properties, derived by 
combining the structural properties of the 156 surviv- 
ing satellites from all eleven halos. The histograms are 
visually similar. (Note that the LMC and SMG are not 
included in the observational data set since they are rota- 
tionally supported and our models are restricted to hot 
systems. They would be equivalent to the most lumi- 
nous, highest velocity dispersion systems in our model 
data set that appear to be missing from the Milky Way 
distribution.) 

To quantify the level of similarity of the simulated and 
observed data s ets we use the 3-dimensional KS-statistic 
(jPasano fc Franccsc hini 1987) 

Z„^3B = dmaxVn, (18) 

where n is the number in the sample tested against our 
model parent distribution of all 156 surviving satellites. 
In this method dmax is defined as the maximum differ- 
ence between the observed and predicted normalized in- 
tegral distributions, cumulated within the eight volumes 
of the three-dimensional space defined for each data point 
{Xi, Yi, Zi) = {^io.i,L-k, i, a^j) by 

(x < Xi,y < Yi,z < Zi), {x > Xi,y > Yi,z > Zi). 

(19) 



iFasano fc Franceschinil l)1987fl present assessments of the 
significance level of values obtained for Z„^3d as a func- 
tion of n and of the degree of correlation of the data. 
Since we already have eleven similarly-sized samples 
drawn from the same parent distribution, we instead 
quantify the significance level of Zn^su found for the 
Milky Way satellites by comparing it to the distribu- 
tion of Zn,3T> for our simulated samples. Figure |H| shows 
a histogram of the results for our simulated halos, with 
the dotted line indicating where the Milky Way satel- 
lite distribution falls. According to this test only one 
of the eleven simulated populations is more similar to 
the simulated parent population than the observed satel- 
lites. (Note that ~80% of our simulated samples have 
Zn,3D < 1.2. Th is significance level is similar to those de- 
rived bv lFasano fc Franccschini (1987) for 3-dimcnsional 
samples with n = 10 and a moderate degree of correla- 
tion in the distribution — see their Figure 7 — as might 
be expected given the expected relation between cto and 
Ltot, see §3.2.2.) 

3.2.2. Correlations in satellite structural parameters 

Figure shows the relationship between the central 
(< Tc), 1-D light-weighted velocity dispersion and satel- 
lite stellar mass, Mi,, for model galaxies and observed 
galaxies in the Local Group. Crosses show surviving 
model satellites for all halos and open circles show the 
relationship for the same set of satellites before they were 
accreted into the host dark matter halo. Solid triangles 
show the relationship for Local Group satellites as com- 
piled by D ekeLfcJVoQ (2003). The two nearly identi- 
cal solid lines show the best-fit regressions for the initial 
and final model populations. The dashed line shows the 
best-fit line for the data. Our model galaxies reproduce 
a trend quite similar to that seen in the data. The rel- 
ative agreement is significant for two reasons. First, the 
stellar velocity dispersion of our initialized satellites is 
set by the underlying potential well of their dark matter 
halos convolved with their associated King profile pa- 
rameters. While in H2.4I we set King profile parameters 
using a phenomenological relation based on the stellar 
luminosity (Li, => fc), there was no guarantee that the 
dark matter potential associated with a given luminos- 
ity would provide a consistent stellar velocity dispersion 
(rc+PoM ^ "'*)■ III this sense, the general agreement be- 
tween model satellites and the data is a success of our star 
formation prescription, which varies based on the mass 
accretion histories of halos of a given size (and therefore 
density structure). 

A second interesting feature shown in Figure |7| is that 
final surviving satellites obey the same relation as the 
initial satellites. Most of these systems have experienced 
significant dark matter mass loss, but since the star par- 
ticles are more tightly bound, their velocity dispersion 
does not significantly evolve. This point is emphasized 
in Figure IHl where we plot the central (< Tc), 1-D ve- 
locity dispersion for the dark matter in halos, again as a 
function of the satellite galaxy's stellar mass, M^,. As in 
Figured open circles show the relationship for the final, 
surviving satellites, and crosses show the relationship for 
those same satellites before they were accreted. Unlike 
in the case of light-weighted dispersions, the dark matter 
dispersion velocities in the surviving systems is system- 
atically lower than in the initial halos owing to the loss of 
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the most energetic particles. They also exhibit a broader 
scatter at fixed stellar mass, reflecting variations in their 
mass- loss histories. Comparing again to Figure [7| we see 
that most of the particles associated with light in these 
systems remains bound to the satellites and their velocity 
dispersions do not evolve significantly. This result may 
have important implications for interpreting the nature 
of the dark matter halos of dwarf galaxies in the Local 
Group and for understanding the regularity in observed 
dwarf properties irrespective of their environments. In 
future work we will return to a more detailed structural 
and evolutionary analysis of the light matter and the 
dark matter halos in which the stars are embedded. 

The results presented in this and the previous sub- 
section clearly indicate that our star formation scenario 
coupled with setting the King parameters of our infalling 
dwarfs to match Local Group observations lead to surviv- 
ing satellite populations consistent both in number and 
structural properties with the Milky Way's. 

3.2.3. The stellar halo 's mass and density profile 

Estimates for the size, shape and extent of the 
Milky Way's stell ar halo come either from star count 
surveys fM orrison "eraI] l2000t IChiba fc Beer j I2OO0I 
Yannv et a l. 2000; Sicgcl. Maie wski. Reid. fc ThompsonI 
or studies where distances could be estimated us- 
ing R R Lyraes l|Wetterer fc McGrawiri996t llvezic et al.1 
I2OOOI) . These studies agree on a total luminosity of or- 
der Ly ~ 10^ Lq (or mass ~ 2 x lO^Af©), which is in 
good agreement with the unbound stellar luminosity for 
all eleven of our model stellar halos, listed in Column 
6 of Table 1 (numbers in brackets again refer to stars 
from accretion events since the last >10% merger). The 
match between predicted and observed total halo mass is 
non-trivial and depends sensitively on the mass accretion 
history of the dark matter halo along with the value of 
the star-formation timescale, t*. Specifically, we show in 
§4.1.1 that the majority of dwarf galaxies that make up 
the stellar halo were accreted early, more than ~ 8Gyr 
ago. The total stellar halo mass (~ lO^M©) is relatively 
small compared to the total cold baryonic mass in ac- 
creted satellites (~ IO^^Mq), because the star formation 
timescale is long compared to the age of the Universe at 
typical accretion times, and the stellar mass fractions are 
correspondingly low (see Figure^J ■ If we would have cho- 
sen a star formation timescale short compared the time 
of typical accretion for a destroyed system (e.g. ^ 5Gyr) 
this would have resulted in a stellar halo of stripped stel- 
lar much more massive than that observed f or the Milky 
Way. T his is in agreement with the results ofBroo k et al.l 
lj2004b(l . who found that a strong feedback model (effec- 
tively slowing the star formation rate in dwarf galaxies) 
in their smoothed particle hydrodynamical simulations 
of galaxy formation was necessary in order to build rela- 
tively small halo components in their models. 

The observational studies find density profiles falling 
more steeply than the dark matter halo (a power law in- 
dex in the range -2.5 to -3.5, compared to ^ —2 for the 
dark matter at relevant scales). Some of the variance 
between results from different groups can be attributed 
to substructure in the halo since these studies have com- 
monly been limited in sky-coverage with surveys covering 
significant portions of the sky only now becoming feasi- 
ble. Figure El plots the density profiles generated (ar- 



bitrarily normalized) from our four representative stellar 
halo models (light solid curves), which transition between 
slopes of -2 within ~ lOkpc to -4 around 50kpc and fall off 
even more steeply beyond this. To illustrate the general 
agreement with observations, the dotted line is a power 
law with exponent of -3. Note that there is some vari- 
ation in the total luminosity (about a factor of 2) and 
slopes of our model halos, as might be expected given 
their different accretion histories. There is also a clear 
roll-over below the power law in the outer parts of the 
stellar halo, sometimes at radii as small as 30kpc. 

To contrast to the light, the density profiles of the dark 
matter in our models are plotted in bold lines Figure El 
(also with arbitrary normalization). The dark matter 
profiles are all close to an NFW profile with mhaio = 
1.4 X lO^^Af© and rhaio = lOkpc. Within - 30kpc of the 
Galactic center it appears that our stellar halos roughly 
track the dark matter, but beyond this they tend to fall 
more steeply. The difference in profile shapes — and 
the steep roll-over in the light matter at moderate to 
large radii — is a natural consequence of embedding the 
light matter deep within the dark matter satellites: The 
satellites' orbits can decay significantly before any of the 
more tightly bound material is lost. Hence we anticipate 
that more/less extended stellar satellites would result in 
a more/less extended stellar halo. Studies of the distant 
Milky Way halo are still sufficiently limited that it is not 
possible to say whether the location of the roll-over in 
our model stellar halos is in agreement with observations 
and this could be an interesting test of our models in the 
near future (see, e.g. Ivczic et al. 2003). 

4. RESULTS II: MODEL PREDICTIONS 

We have now fixed our free parameters to be Zro = 10, 
t-i, = 15Gyr and /gas=0.02. By limiting our description 
of the evolution of the baryons associated with each dark 
matter satellite to depend on only these parameters, we 
find we have little freedom in how we choose them. For 
example, if we were to choose a shorter star formation 
timescale i*, we would over-produce the mass of the stel- 
lar halo, form dwarf galaxies that were over-luminous 
at fixed velocity dispersions, and form dwarfs with low 
gas fractions compared to isolated dwarfs observed in the 
Local Group. The first two problems could be adjusted 
by adapting /gas, but the last problem is independent of 
this. 

Despite its simplicity, our model reproduces observa- 
tions of the Milky Way in some detail. In particular, 
we recover the full distribution of satellites in structural 
properties. This suggests both that have we assigned the 
right fraction of dark matter halos to be luminous and 
that our luminous satellites are sitting inside the right 
mass dark matter halos. 

We can now go on with some confidence to discuss the 
implications of our model for the mass accretion history 
of the halo and satellite systems (§4.1) and the level of 
substructure in the stellar halo (§4.2). 

4.1. Building up the stellar halo and satellite systems 

4.1.1. Accretion times and mass contributions of 
infalling satellites 

The stellar halo in our model is formed from stars orig- 
inally born in accreted satellites. Once accreted, satel- 
lites lose mass with time until the satellite is destroyed. 
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Once a particle becomes unbound from a satellite, we 
associate its stellar mass with the stellar halo. Figure 
[TUI shows the cumulative luminosity fraction of the stel- 
lar halo (solid lines) coming from accreted satellites as 
a function of the accretion time of the satellite for halos 
1,2, 6, and 9. Clearly most of the mass in the stellar ha- 
los originates in satellites that were accreted more than 
8 Gyr ago. The dotted lines show the contribution to 
the stellar halo from satellite halos more massive than 
Mvir > 2 X IO^^Mq at the time of their accretion. While 
only 10 — 20 of the ~ 150 accreted satellites meet this 
mass requirement, we see that ~ 75 — 90% of the mass as- 
sociated with each stellar halo originated within massive 
satellites of this type. 

Compare these to the dashed lines, which show the cu- 
mulative number fraction of surviving satellite galaxies 
as a function of the time they were accreted for the entire 
population (long-dashed lines) and restricted to satellite 
halos that were more massive than Mvir^ 5 x IO^Mq at 
the time of their accretion (short dashed lines). We see 
that surviving satellites are accreted much later (~ 3 — 5 
Gyr lookback) than their destroyed counterparts and 
that the most massive satellites that survive tend to be 
accreted even later because the destructive effects of dy- 
namical friction are more important for massive satel- 
lites. 

4.1.2. Spatial growth 

Studies of dark matter halos in N-body simulations 
show that they are built fro m the inside out (e.g. 
IHelmi. White, fc Springeill2003|) . The top panel of Fig- 
ure ^2 confirms that this idea holds for our model stel- 
lar halos: it shows the average over all our halos of the 
fraction of material in each spherical shell from all accre- 
tion events (solid line) and from those that have occurred 
since the last major (> 10%) merger (dotted line — the 
time when this occurred is given in column 4 of Table 
nj. Although the recent events represent only a fraction 
of the total halo luminosity (~5-50%, see Table they 
become the dominant contributor at radii of 30-60kpc 
and beyond. There is some suggestion of this being the 
case for the Milky Way's halo globular cluster popula- 
tion, which can can fairly clearly be separated into an 
'old' , inner population (which exhibits some rotation, 
is slightly flattened and has a metallicity gradient) and 
a 'young' outer one (which is more extended and has a 
higher velocity dispersion — see Zinn 1993). 

One implication of the inside-out growth of stellar ha- 
los, combined with the late accretion time of surviving 
satellites is that the two should each follow different ra- 
dial distributions. The dashed lines in the bottom panel 
of Figure shows the number fraction of all surviving 
satellites in our models as a function of radius — the 
distribution is much flatter than the one shown for the 
halo in the upper panel. In fact, all satellites of our own 
Milky Way (except Sgr) lie at or beyond 50kpc from its 
center, with most 50-150kpc away, — as shown by the 
solid line in the lower panel. Hence, the radial distri- 
bution of the observed satellites is consistent with our 
models and suggests that they do indeed represent re- 
cent accretion events. 

4.1.3. Implications for the abundance distributions of 
the stellar halo and satellites 



Studies which contrast abundance patterns and stel- 
lar populations in the stellar halo with those in dwarf 
galaxies seem to be at odds with models (such as 
ours) that build the stellar halo from satellite accre- 
tion (Una vane. Wvse. fc Gilmorc"1996'). For example, 
both field and satellite populations have similar metal- 
licity ranges, but the former typicall y have higher alpha- 
element abundances than the latter llTolstov et alJl2003t 
IShetrone et a,l 1120031 1 Venn et 31.1120041). Clearlv. it is not 
possible to build the halo from present-day satellites. 

We have already shown (in Figure ^1 ) that we would 
expect a random sample of halo stars to have been 
accreted 8-lOGyears ago from satellites with masses 
Mvir^ 1O^°M0, while surviving satellites are accreted 
much more recently. (Note that Figure E| deliberately 
compares the cumulative luminosity fraction of the stellar 
halo to the 7T,Mm6er fraction of satellites. This is the most 
relevant comparison to make when interpreting observa- 
tions because any sample of halo stars will be weighted by 
the luminosity of the contributing satellites, while sam- 
ples of satellite stars are often composed of a few stars 
from each satellite.) Figure IT^ explores the number and 
luminosity contribution of different luminosity satellites 
to each population in more detail. It shows the number 
fraction of satellites in different luminosity ranges con- 
tributing to the stellar halo (dotted lines) and satellites 
system (dashed lines): the peak of the dotted/dashed 
lines at lower/higher luminosities is a reflection of the 
much later accretion time — and hence longer time avail- 
able for growth of the individual contributors — of the 
satellite system relative to the stellar halo. However, as 
discussed above, it is more meaningful to compare the 
number fraction of surviving satellites to the luminosity 
fraction of the halo (solid lines) contributed by satellites 
of a given luminosity range. The solid line emphasizes 
(as noted above) that most of the stellar halo comes from 
the few most massive (and hence most luminous) satel- 
lites, with luminosities in the range 10^ — lO^L©. In 
contrast. Galactic satellite stellar samples would likely 
be dominated by stars born in 10^ — lO^i© systems. 

Overall our results provide a simple explanation of the 
difference between halo and satellite stellar populations 
and abundance patterns. The bulk of the stellar halo 
comes from massive satellites that were accreted early, 
and hence had star formation histories that must be short 
(because of their early disruption) and intense (in order 
to build a significant luminosity in the time before dis- 
ruption). In contrast, surviving satellites are lower mass 
and accreted much later, and hence have more extended, 
lower level star formation histories. Stars formed in these 
latter environments represent a negligible fraction of the 
stellar halo in all our models. This is confirmed by the 
last column of Table which lists the percentage contri- 
bution of surviving satellites to the total halo (less than 
10% in every case). Note that the contributions of sur- 
viving satellites to the local halo (i.e. within 10-20kpc of 
the Sun) , which is the only region of the halo where de- 
tailed abundance studies have been performed, are even 
lower (less than 1% in every case). 

A more quantitative investigation of the consequences 
of the difference between the " accr etion age" of stars and 
satellites in the h alo is underway ijRobertson et al.ll2005t 
iFont et alJl2005|) . 
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4.2. Substructure 

Abundant substructure is one of the most basic expec- 
tations for a hierarchically-formed stellar halo. Here we 
give a short description of the substructure we see in our 
simulations, and reserve more detailed and quantitative 
explorations for future work. Recall, our study (by de- 
sign) follows the more recent accretion events in our halo 
more accurately than the earlier ones — we showed in 
§4.1.2 that these are the dominant contributors to the 
halo at radii of 30-60kpc and beyond. Hence we can ex- 
pect our study make fairly accurate predictions of expec- 
tations of the level of substructure in the outer parts of 
galactic halos — precisely the region where substructure 
should be more dominant and easier to detect. 

Figures El and El show external galactic views for halo 
realizations 1, 2, 6, and 9. The color code reflects surface 
brightness per pixel: white, 24 magnitudes per square 
arcsecond, to light blue at 30 magnitudes per square arc- 
second, to black which is (fainter than) 38 magnitudes 
per square arcsecond. The darkest blue features are of 
course too faint to be seen (except by star counts). We 
have simply set the scale in order to reveal all the spatial 
features that are there in principle. We mention that our 
test particles f H2.2.4|l were used in making these images. 

In addition to spatial structure, we also expect signifi- 
cant structure in phase space. A two-dimensional slice of 
the full six dimensional phase space is illustrated in Fig- 
ure El where we plot radial velocity Vr versus radius r 
for all of halo 1 (left) and halo 9 (right). Each point rep- 
resents 1000 solar luminosities. The color code reflects 
the time the particle became unbound from its original 
satellite: dark blue for particles that became unbound 
more than 12 Gyr ago and white for particles that either 
remain bound or became unbound less than 1.5 Gyr ago. 
The radial gradient in color reflects the "inside out" for- 
mation of the stellar halo discussed in previous sections. 

Note that significant coherent structure is visible in 
Figure [TBI even without any spatial slicing of the halos. 
Except for particles belonging to bound satellites (white 
streaks), the structure strongly resembles a nested se- 
ries of orbit diagrams. This is not surprising since the 
halo was formed by particles brought in on satellite or- 
bits. A direct test of this prediction should be possible 
with SDSS data and other similar surveys. Indeed, if 
the phase space structure of stellar halo stars does reveal 
this kind of orbit-type structure it will be a direct indi- 
cation that the stellar halo formed hierarchically. Figure 
1161 shows the same diagram for halo 1 , now subdivided 
into two distinct quarters of the sky. 

5. SUMMARY AND CONCLUSIONS 

We have presented a cosmologically self-consistent 
model for the formation of the stellar halo in Milky- 
Way type galaxies. Our approach is hybrid. We use a 
semi-analytic formalism to calculate a statistical ensem- 
ble of accretion histories for Milky Way size halos and to 
model star formation in each accreted system. We use a 
self-consistent N-body approach to follow the dynamical 

^ In most cases, we subsample our luminous particles in order to 
plot a single point for every IOOOLq . However, some of the particles 
in our simulations have luminosity weights greater than 1000 Lq. 
In these cases we plot a number of points (= iparticlc/lOOO) with 
the same Vr and r as the relevant particle (using small random 
offsets to give the effect of a "bigger" point on the plot). 



evolution of the accreted satellite galaxies. A crucial in- 
gredient in our model is the explicit distinction between 
the evolution of light and dark matter in accreted galax- 
ies. Stellar material is much more tightly bound than 
the majority of the dark matter in accreted halos and 
this plays an important role in the final density distribu- 
tion of stripped stellar material as well as the evolution 
in the observable quantities of satellite galaxies. 

A primary goal of this, our first paper on stellar halo 
formation, was to normalize our model to, and demon- 
strate consistency with, the gross properties of the Milky 
Way stellar halo and its satellite galaxy population. We 
constrained our two main star formation parameters, 
the redshift of reionization Zro, and the star formation 
timescale of cold gas t^, using the observed number 
counts of Milky Way satellites and the gas mass fraction 
of isolated dwarf galaxies. With these parameters fixed, 
the model reproduces many of the observed structural 
properties of (surviving) Milky Way satellites: the lumi- 
nosity function, the luminosity-velocity dispersion rela- 
tion, and the surface brightness distribution. The satel- 
lite galaxies that are accreted and destroyed in our model 
produce stellar halos of material with a total luminosity 
in line with estimates for the stellar halo of the Milky 
Way (~ lO^io). 

The success of our model lends support to the hierar- 
chical stellar halo formation scenario, where the stellar 
halos of large galaxies form mainly via the accretion sub- 
sequent disruption of smaller galaxies. More specifically, 
it allows us to make more confident predictions concern- 
ing the precise nature of stellar halos and their associated 
satellite systems in Milky Way type galaxies. These in- 
clude: 

• The density profile of the stellar halo should fol- 
low a varying power-law distribution, changing in 
radial slope from ~ —2 within ~ 20kpc to ~ —4 
beyond 50kpc. The distribution is expected to be 
much more centrally concentrated than the dark 
matter, owing to the fact that the stars that build 
the stellar halo were much more tightly bound to 
their host systems than the dark material respon- 
sible for building up the dark matter halo. 

• Stellar halos (like dark matter halos) are expected 
to form from the inside out, with the majority of 
mass being deposited from the ~ 15 most massive 
accretion events, typically dwarf-irregular size ha- 
los with mass ~ 10^°Mq and luminosities of order 

lO^-lO^ig. 

• Destroyed satellites contributing mass to the stellar 
halo tend to be accreted earlier than satellites that 
survive as present-day dwarf satellites (~ 9Gyr 
compared to ^ 5 Gyr in the past). 

• Substructure, visible both spatially and in phase 
space diagrams, should be abundant in the outer 
parts of galaxies. Proper counts of this structure, 
both in our galaxy and external systems, should 
provide important constraints on the late-time ac- 
cretion histories of galaxies and a test of hierarchi- 
cal structure formation. 

Together, the second and third points imply that most 
of the stars in the inner halo are associated with massive 
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satellites that were accreted 9Gyr ago. Dwarf satel- 
lites, on the other hand, tend to be lower mass and are 
associated with later time accretion events. This sug- 
gests that classic "stellar halo" stars should be quite dis- 
tinct chemically from stars in surviving dwarf satellites. 
We explore this point further in two companion papers 



(Robertson et al. 05; Font et al. 05). 
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TABLE 1 

Properties of our simulated stellar halos. 

Note. — Numbers in brackets in columns 5, 6 and 7 are for events since the last > 10% merger (the time of this 
event is given in column 4). 
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Fig. 1. — The mass assembly histories for halos 1, 2, 6, and 9. SoUd jagged lines show the histories for the dark matter halos, Mvir{< t), 
generated via the EPS merger tree method. Smooth solid lines show our best-fit "smooth" accretion history used for the background 
potential in our N-body simulations (Equation |^. Smooth dotted and short-dashed lines show the evolution in our disk and bulge 
component masses used to set the galactic contributions to the N-body potentials (Equations The long-dashed lines show our (main) 

results for stellar halo assembly histories. Stellar mass is assigned to the stellar halo component at the time it becomes unbound from an 
accreted satellite halo. Note that the stellar halo in realization 6 is relatively small compared to the other systems whose dark matter halos 
are all the same size. This is because dark matter halo 6 accumulated most of its mass relatively early, when the stellar mass fractions of 
accreted satellites were small. Note also that halo 9 has experienced a very recent, massive disruption accretion event, which caus es its 
stellar halo mass to increase sharply at a recent lookback time. This event is seen clearly in the image of halo 9 shown in Figure [T^ 
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Fig. 2. — The energy distribution function of our initial condition dark matter halo (histogram, dM/de) along with three example energy 
distributions for stellar matter, (dAf/de)*, in satellites. The mass to light ratio of each particle of energy e is assigned based on the ratio 



of {dM/de)i, to (dM/de). Energy in this plot is in units of GAf|g/2Rijaio- 
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Fig. 3. — V-band luminosity plotted against the ratio of mass in HI to V-band luminos ity for satellites of the Milky Way and Andromeda 
(open symbols) and field dwarfs (solid symbols). Data is taken from the compilation bv lGrebel. Gallagher. Ik Harbec^ 120031) 




Fig. 4. — Plotted is the gas mass per stellar luminosity ratio Mgaa/iv for each accreted satellite at the time of their accretion, taccn for 
four example halos. Early accreted systems are more gas rich owing to the rather long star formation timescales in these systems. Solid 
points points indicate satellites that survive until the present day. 
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Fig. 5. — Histograms of the distribution of Milky Way satellites (solid lines) and simulated satellites (dashed lines — from all 11 halos) 
as functions of observed quantities 
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Fig. 6. — Derived KS statistic for the satellite distribution in each simulated halo compared with those from the combined sample (solid 
line) and the observed Milky Way distribution (dotted line). 
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Fig. 7. — The relationship between the central, 1-D stellar velocity dispersion and satellite mass for all of our surviving model satellites. 
Crosses show surviving satellites at the current epoch and open circles show the relationship computed before the satellites were accreted 
into the host dark matter halo. Solid triangles show the same relationship for local group satellites as compiled by Dckcl & Woo (200.3). 
The two nearly identical solid lines show the best-fit regressions for the initial and final model populations. The dashed line shows the 
best-fit line for the data. The best-fit lines all have similar slopes, cr, oc A/" with a ~ 0.2 (a = 0.19, 0.18 and 0.19 for initial, final, and data 
populations respectively, and with errors of 0.01 in a and ~ 0.07 in logarithmic normalization.). The model satellites match the observed 
trend in the data quite well, considering the observational uncertainties (see, e.g. the discussion in lDekel fc Wooll2003l) . 




Fig. 8. — The the central, 1-D velocity dispersion in the dark matter versus satellite stellar mass. As in Figure|7| crosses show surviving 
model satellites and open circles show the relationship computed before those systems were accreted into the host dark matter halo. The 
lines again show least square regression fits, ctdm oc MJ. The final systems tend to have lower velocity dispersions and a broader scatter 
at fixed stellar mass owing to the loss of loosely-bound dark matter particles after accretion. Compare this result to Figure Q where very 
little shift occurs in the more tightly bound stellar material. 
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Fig. 10. — The cumulative fraction of stellar halo mass built from accreted satellite galaxies as a function of the accretion time of 
the satellites for halos 1-4. Solid lines show the full halo and dotted lines show the contribution from satellite halos more massive than 
M^ir > 2 X 1Q^°Mq at the time of their accretion. While for halos (1,2,6, and 9) only (18,10,12,13) of the (115,102, 182, 153) accreted 
luminous satellites were more massive than 2 X IO^^Mq, we see that ~ 75 — 90% of the mass associated with each stellar halo originated 
within massive satellites of this type. For comparison, the dashed lines show the cumulative fraction of surviving satellite galaxies as a 
function of the time they were accreted. Short-dashed lines also show the cumulative accretion times of surviving satellites, except now 
restricted to satellite halos that were more massive than Mvir^ 5 X IO^Mq at the time of their accretion We see that surviving halos tend 
to be accreted later than destroyed halos. There is also a tendency for massive satellites that survive to be accreted even later because the 
destructive effects of dynamical friction are more important for massive satellites. 
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Fic;. 11. — Top panel compares the luminosity fraction of the entire stellar halo spherical shells at radius r from the center of the parent 
galaxy (solid lines) with the contribution from events since the last > 10% merger event (dotted lines). Bottom panel compares the number 
fraction of Milky Way satellites in spherical shells as a function of Galactocentric distance (solid line) with the number fraction of our 
surviving satellite population taken from all eleven simulated halos (dashed line). 
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Fig. 12. — The number fraction of all events, binned In satellite lunminosity, contributing to all eleven of our simulated halos (dotted 
lines) compared to the number fraction of surviving satellites (dashed lines). The solid lines show the fraction of the total luminosity 
contributed by each range of satellite luminosities. 
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Fig. 13. — "External Galaxy" views for halo realizations 1 (left) and 2 (right). The boxes are 300 kpc by 300 kpc. The blue/white color 
scale indicates surface brightness: 23 Magnitudes per square arcsecond (white) to 38 Magnitudes per square arcsecond (dark blue / black) 
where we have assumed a stellar mass to light ratio of 2. The eye picks up lighter blue (middle of the bar) at about 30 Magnitudes per 
square arcsecond. 




Fig. 14. — "External Galaxy" images for halos 6 (left) and 9 (right). The color codes are the same as those in Figure ITsl A recent 
disruption has occurred in halo 9 (~ l.SGyr lookback time) and the residue of this event is seen as the bright plume running from the 
"north-west" of the halo (upper left) down towards the halo center. The bright feature just to the "south-west" of halo 9's center is also 
associated with the same disruption event. 
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Fig. 15. — Radial phase space diagrams {Vr vs r relative to the host halo center) for halos 1 (left) and 9 (right). Each point represents 
1000 solar luminosities. The color code reflects the time each particle became unbound to its parent satellite. White points are either 
bound or became unbound in the last 1.5 Gyr, while dark blue points became unbound more than 12 Gyr ago. The radial color gradient 
reflects the tendency for inner halo stars to be accreted (and stripped) early in the Galaxy's history. The white feature at r ~ 80 kpc in 
halo 9 represents a very recent disruption event - the most recent massive disruption seen in our ensemble of 11 halo realizations. 




Fig. 16. — Radial phase space diagrams for halo 2, w here the left and right figures represent to separate quarters of the sky. The color 
code and a,xis labels are the same as those in Figure [T^ 



